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A  BARELY  IMPLICIT  CORRECTION  FOR  FLUX-CORRECTED  TRANSPORT 

L  INTRODUCTION 

The  solution  of  time  dependent  compressible  flow  problems  is  complicated 
by  conflicting  requirements  of  mathematical  accuracy,  nonlinearity,  physical 
conservation  and  positivity.  This  is  especially  true  near  discontinuities  where 
"accurate”  high  order  algorithms  produce  ripples  while  linear  monotonic  (i.e. 
positivity  preserving)  schemes  are  highly  diffusive.  After  Godunov  [1]  showed 
that  a  linear  algorithm  ensures  positivity  only  if  it  is  first  order,  the  next  logical 
step  was  to  look  at  nonlinear  methods  to  develop  effectively  higher  order,  more 
accurate  monotonic  schemes.  The  first  high-order  monotone  algorithm  (Boris 
[2])  was  designed  to  maintain  local  positivity  near  steep  gradients  while  keeping 
a  high  order  of  accuracy  elsewhere.  The  major  principles  of  the  monotone  high 
order  algorithms  are  that  they  maintain  positivity  through  a  procedure  that 
uses  a  nonlinear  combination  of  diffusive  and  antidiffusive  fluxes.  The  Flux- 
Corrected  Transport  (FCT)  algorithm  that  we  use  in  this  paper  [3,  4]  is  made 
fourth  order  by  the  appropriate  subtraction  of  corrected  fluxes.  Other  mono¬ 
tone  methods  have  been  reviewed  by  Woodward  and  Collela  [5]  and  Baer  [6]. 

In  this  paper  we  confine  our  discussions  to  a  Barely  Implicit  Correction  (BIC) 
to  FCT.  BIC  is  also  extendable  to  other  monotone  methods. 

Positivity-preserving  monotone  FCT  methods  were  developed  to  calculate 
shocks  accurately.  Even  for  subsonic  flows  with  discontinuities,  their  high  accu¬ 
racy  produced  much  better  solutions  than  standard  finite  difference  techniques. 

The  early  FCT  methods  were  explicit.  No  serious  limitation  arose  from  the 
explicitness  in  supersonic  flows  because  the  major  features  of  interest  in  the 
flow  move  at  about  the  sound  speed.  Using  these  methods  for  subsonic  flows, 
however,  is  economical  only  if  the  characteristic  velocities  in  the  flow  field  are 
a  reasonable  fraction  of  the  speed  of  sound  [7,  8]  or  if  the  fast  sound  waves  are 
mathematically  removed  from  the  system  of  equations. 
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The  Barely  Implicit  Correction  described  in  this  paper  was  motivated  by 
the  need  to  calculate  subsonic  flows  accurately  in  which  the  velocities  of  the 
important  flow  structures  are  much  lower  than  the  speed  of  sound.  In  typical 
cases,  we  are  interested  in  flow  velocities  from  centimeters  to  tens  of  meters 
per  second.  These  flow  velocities  are  encountered,  for  example,  in  laminar 
flames  and  low-speed  fuel  injection  in  engines.  Our  objectives  are  to  remove  the 
timestep  limit  imposed  by  the  speed  of  sound,  retain  the  accuracy  required  to 
resolve  the  detailed  features  of  the  flow,  and  reduce  the  computational  costs. 

The  obvious  way  to  beat  the  sound-speed  limit  on  the  timestep  is  to  make 
the  calculation  implicit.  This  has  been  done  successfully  for  many  linear  meth¬ 
ods,  such  as  the  MacCormack  method  [9],  the  Beam  and  Warming  method  [10], 
and  the  ICE  and  RICE  methods  of  Hirt  and  Cook  [11].  In  addition,  recent  de¬ 
velopments  have  been  reported  for  implicit,  nonlinear  PPM  [12]  and  TVD  [13] 
methods.  A  major  problem  with  all  of  these  methods  is  that  they  are  relatively 
expensive,  even  though  they  can  be  made  relatively  accurate. 

Another  approach  is  the  asymptotic  methods.  Examples  of  these  are  the 
methods  developed  by  Jones  and  Boris  [14],  Rehm  and  Baum  [15]  and  Paolucci 
[16].  In  these  methods,  the  only  effects  of  compression  that  are  allowed  are  the 
changes  in  density  due  to  heating  or  cooling.  Pressure  fluctuations  are  filtered 
out,  thus  removing  the  timestep  limit  imposed  by  the  sound  speed.  However, 
other  effects  from  sound  waves  are  removed  in  this  process. 

A  a  useful  approach  was  given  by  Casulli  and  Greenspan  [17].  Their  analysis 
indicated  that  it  is  not  necessary  to  treat  all  of  the  terms  in  the  gas  dynamic 
equations  implicitly  to  be  able  to  use  longer  timesteps  than  those  dictated  by 
explicit  stability  limits.  Only  those  explicit  terms  which  force  this  limit  need  to 
be  treated  implicitly. 

The  algorithm  presented  in  this  paper  has  two  steps.  The  first  step  is 
explicit.  It  is  performed  at  a  large  timestep  governed  by  a  CFL  condition  on 
the  fluid  velocity.  This  step  should  be  done  with  an  accurate  nonlinear  monotone 
method,  and  we  have  used  FCT  in  the  examples  given.  The  second  step  is  an 
implicit  correction  step  requiring  the  solution  of  one  elliptic  equation  for  the 
pressure  correction.  The  term  “Barely  Implicit  Correction”  emphasizes  our  use 


of  the  idea  of  Casulii  and  Greenspan,  that  only  certain  terms  must  be  treated 
implicitly. 

The  total  cost  per  timestep  of  BIC-FCT  is  about  the  same  as  for  a  full 
explicit  FCT  step.  Thus  the  cost  of  a  complete  calculation  is  one  or  two  orders 
of  magnitude  below  that  required  if  a  very  slow  flow  were  treated  explicitly. 
Since  only  one  elliptic  equation  is  solved,  the  method  is  considerably  faster 
than  many  implicit  methods  commonly  used.  In  addition,  using  a  nonlinear 
monotone  method  for  the  explicit  step  ensures  high  accuracy. 

fl.  METHOD  OF  SOLUTION 

Derivation  of  the  Barely  Implicit  Correction 

We  are  solving  the  compressible  gas  dynamics  conservation  equations  for 
density  p,  momentum  density  pv,  and  total  internal  energy,  E, 

^  =  -V/xv,  (1) 

=  -v-rv  V-VP,  (2) 

=  -V  •  (B  +  P)v  ,  (3) 

where  the  total  energy  density  E  is 

E  =  e+^pv2.  (4) 

The  equation  of  state  relating  pressure  and  internal  energy  is 

P  =  (7  -  l)e  •  (5) 

In  a  recent  paper,  Casulii  and  Greenspan  [17]  showed  that  it  is  not  necessary  to 
treat  every  term  in  a  finite-difference  algorithm  implicitly  to  avoid  the  timestep 
constraint  imposed  by  the  Courant  condition.  Further,  they  showed  that  only 
those  terms  containing  the  pressure  in  Eq.  (2)  and  the  velocity  in  Eq.  (3)  must 
be  treated  implicitly.  Their  paper  provides  the  starting  concepts  for  the  work 
we  present.  In  addition,  we  have  extended  their  analysis  to  include  an  implic- 


itness  parameter,  w,  that  can  be  used  to  vary  the  degree  of  Implicitness  of  the 
algorithm.  In  general,  we  can  have  0.5  <  w  <  1,  where  the  implicit  terms  are 
centered  in  time  for  oj  =  0.5.  For  uj  <  0.5,  the  method  is  found  to  be  unstable 
for  sufficiently  large  timesteps. 

There  are  two  stages  to  the  algorithm.  One  stage  is  an  explicit  predictor 
that  determines  the  provisional  values  p  and  v, 


P-P 


pv  -  p°v° 


=  -  V  •  p°y°  , 

=  —  V  •  p°v0y°  -  VP®  , 


The  tilde  denotes  predictor  values  at  the  new  time,  and  the  superscripts  o  and 
n  are  used  to  denote  the  old  time  and  new  time,  respectively.  So  far  only  time 
has  been  differenced,  not  space.  The  implicit  forms  of  Eqs.  (2)  and  (3)  are 


pn\n  —  p°\° 

At 

En  -  E° 


=  -V  ■  p°w°y°  -  V[o/Pn  +  (1  -  w)P“l  , 
=  -V  •  ( E°  +  P°)\u)vn  +  (1  -  u/)v°]  , 


where  w  is  the  implicitness  parameter  discussed  above.  When  u  =  1,  the  algo¬ 
rithm  is  completely  implicit  and  reverts  to  the  original  equations  analyzed  by 
Casulli  and  Greenspan. 

We  can  reduce  this  implicit  system  to  only  one  equation  by  eliminating  vn 
between  Eq.  (8)  and  (9).  To  do  this,  we  first  define  the  change  in  pressure,  6Py 
as 

6P  =  uj{Pn-P°)  .  (10) 

Then  the  correction  equation  for  momentum  can  be  obtained  in  terms  of  6P  by 
subtracting  Eq.  (7)  from  Eq.  (8), 


pn\n  -  pv 


=  -Vu/(Pn  -  P°)  =  -V6P. 


We  obtain  the  new  velocity  by  rearranging  Eq.  (ll)  and  letting  pn  =  p ,  so 


v"  =  -~V6P  +  v  . 
P 


1 


We  obtain  a  correction  equation  for  energy  using  the  equation  of  state  with  7 
constant,  Eq.  (4), 


£n  = 


SP 


(7  -  l)o  +  e°  ’  (13) 

where  the  oj  factor  appears  from  the  definition  of  6P.  We  find  6P  by  substituting 
Eqs.  (12)  and  (13)  into  Eq.  (9), 


p\2  -  p°\°2 

2A t 


+ 


6P 


=  wAt  v  •  vsp 


(7  -  l)wAt  \  p 

-  wV  •  ( E°  +  P°)v 


(14) 


-  (1  —  w)V  •  (E°  +  P°)v°  . 


Note  that  the  kinetic  energy  change  is  included  explicitly.  For  convenience,  we 
define  the  quantity  E, 


E-  E° 
At 


—  V  ■  (E°  +  P°)  [wv  +  (1  —  w)v°  . 


(15) 


This  allows  us  to  rewrite  Eq.  (14), 


6P 

(7  —  l)«At 


—  uiAt  V  • 


E  -  E°  p\2  -  p°\°2 
At  2At 


(16) 


which  provides  us  with  an  elliptic  equation  for  SP.  The  right  hand  side  of 
Eq.  (16)  is  evaluated  explicitly  using  Eq.  (15).  After  the  elliptic  equation  is 
solved  for  SP,  momentum  and  energy  are  corrected  by  Eqs.  (11)  and  (13). 
Note  that  we  started  with  two  equations  with  implicit  terms,  and  now  we  have 
reduced  it  to  one  equation,  Eq.  (16). 

The  Barely  Implicit  Correction  is  carried  out  in  three  stages.  In  the  first, 
Eqs.  (6),  (7),  and  (15)  are  integrated  with  any  one-step  explicit  method.  The 
pressure  correction  equation,  Eq.  (16)  is  solved  by  an  elliptic  solver  in  the  second 
stage.  The  last  stage  requires  the  use  of  Eqs.  (11)  and  (13)  to  obtain  the  final 
values  of  momentum  and  energy  at  the  new  timestep. 


Solution  Procedure 


The  derivation  given  above  does  not  involve  any  specific  choice  of  method 
for  differencing  the  spatial  derivatives.  The  only  restriction  so  far  is  that  the 
spatial  derivatives  must  be  evaluated  at  the  appropriate  time  levels  indicated 


by  the  superscripts.  This  allows  great  flexibility  in  the  choice  of  the  differencing 
scheme  for  these  terms.  Thus  we  can  integrate  the  explicit  predictor  equations, 
Eqs.  (6),  (7),  and  (15)  with  FCT.  This  gives  us  the  benefits  of  using  a  high- 
order  monotone  method.  We  have  given  the  name  BIC-FCT  to  this  particular 
combination  of  BIC  and  FCT.  Tests,  such  as  those  presented  below,  indicate 
that  it  has  the  same  accuracy  and  flexibility  as  FCT. 

At  each  timestep,  the  solution  procedure  we  have  implemented  is  divided 
into  the  three  stages: 

1.  Explicit  predictor  stage 

The  density  and  momentum  are  advanced  explicitly  as  specified  by  Eqs.  (6)  and 
(7)  using  FCT.  This  produces  the  intermediate  quantities,  p  and  pv.  The  v  is 
found  from  pv/p.  Then  v  is  used  to  obtain  E  given  by  Eq.  (15).  FCT  is  also 
used  to  obtain  E. 

2.  Solution  of  Eq.  (16)  for  6P 

In  one  dimension,  the  solution  to  the  difference  form  of  Eq.  (16)  requires  the 
solution  of  a  system  of  linear  equations  by  a  tridiagonal  matrix  solver.  In  two 
dimensions,  the  solution  requires  an  elliptic  solver.  For  the  two-dimensional 
calculations  shown  below,  we  used  a  multigrid  method  [18].  A  substantial  part 
of  the  computer  time  required  in  this  stage  is  in  setting  up  the  coefficients  for 
an  elliptic  equation  solver. 

3.  Momentum  and  energy  corrections 

These  corrections  are  obtained  from  the  pressure  change  6P  using  Eqs.  (11)  and 
(13),  respectively.  These  corrected  values  and  the  density  obtained  explicitly  in 
the  first  stage  are  the  starting  conditions  at  the  new  timestep. 

These  three  stages  are  carried  out  at  every  timestep.  The  derivatives  in  the 
pressure  difference  equation,  Eq.  (16),  are  approximated  by  central  differences. 
All  physical  quantities  are  calculated  at  cell  centers,  and  those  values  needed  at 
cell  interfaces  are  obtained  by  averaging. 

This  technique  can  be  implemented  in  one,  two,  or  three  dimensions.  In 
one  and  two  dimensions,  several  different  geometries  are  possible.  For  example, 
we  have  implemented  two-dimensional  planar  and  axisymmetric  geometries,  and 
one-dimensional  Cartesian,  cylindrical,  and  spherical.  Any  boundary  conditions 


that  are  commonly  used  with  the  standard  FCT  modules  can  be  used  with  this 
algorithm  [7]. 

in.  TESTS  OF  THE  METHOD 

Advection  of  a  One-Dimensional  Contact  Discontinuity 

The  problem  we  consider  first  is  the  flow  of  air  through  a  duct  in  one 
dimension.  The  duct  is  initially  filled  with  air  at  standard  temperature  and 
pressure.  Then  cold  air  with  twice  the  density  flows  into  the  duct.  There  is 
a  contact  discontinuity  at  the  location  where  the  cold,  dense  air  and  normal 
air  meet.  In  the  absence  of  diffusive  processes,  the  contact  discontinuity  should 
move  at  the  velocity  of  the  incoming  air.  This  numerical  test  shows  the  ability 
of  BIC-FCT  to  propagate  a  contact  discontinuity  with  the  same  accuracy  as 
FCT. 

The  computational  domain  was  divided  into  200  evenly  spaced  cells  of  1  cm. 
Initially,  the  discontinuity  was  0.1  m  from  the  inlet.  The  flow  velocity  of  air 
in  the  duct  was  10  m/s.  The  inlet  conditions  corresponding  to  the  cold  air  are 
held  fixed  throughout  the  calculation. 

The  timestep  used  in  this  calculation  is  0.5  ms,  which  is  the  time  required 
for  the  fluid  to  cross  half  a  cell.  This  should  be  compared  to  the  Courant  limit 
of  24  fi s.  Typical  explicit  methods  are  forced  to  employ  a  timestep  of  less  than 
half  of  the  Courant  limit  to  control  the  growth  of  of  perturbations  in  pressure 
and  velocity.  In  this  example,  a  factor  of  forty  to  fifty  is  gained  over  the  explicit 
timestep. 

Figure  1  shows  the  density  profiles  at  intervals  of  fifty  steps.  The  disconti¬ 
nuity,  initially  across  one  cell,  spreads  to  three  or  four  cells  as  it  moves  across  the 
system.  Most  important,  however,  is  that  the  discontinuity  spreads  no  further 
throughout  tne  course  of  the  solution  and  there  are  no  ripples  in  the  solutions. 
Both  of  these  features  are  in  the  underlying  explicit  FCT  algorithm.  The  results 
presented  in  the  figure  were  obtained  with  u>  =  1.  The  influence  of  sound  waves 
in  this  problem  are  negligible,  so  that  any  stable  value  of  w  gives  the  same  result. 


Sound  Wave  Dumping 


In  most  finite-difference  methods,  high-frequency  sound  waves  are  attenu¬ 
ated.  Implicit  methods,  however,  tend  to  damp  all  frequencies,  with  the  lower 
frequencies  damped  least.  The  problem  we  now  present  tests  the  sound-wave 
damping  in  BIC-FCT. 

We  consider  a  closed,  one-dimensional  pipe  1  m  long  in  which  the  fluid 
velocity  was  initialized  with  a  sinusoidal  variation.  The  maximum  amplitude 
of  the  variation  was  1  m/s  at  the  center  of  the  pipe.  Effectively,  the  initial 
conditions  correspond  to  a  sound  wave  in  the  pipe  with  a  wavelength  of  2  m. 

Each  curve  in  Fig.  2  shows  the  fluid  velocity  at  the  center  of  the  pipe  as 
a  function  of  the  number  of  cycles  for  a  different  value  of  u>.  The  damping  is 
greatest  when  u>  =  1,  which  is  when  the  method  is  completely  implicit.  The 
damping  decreases  as  w  is  reduced,  and  it  becomes  negligble  when  tu  =  0.5.  Any 
further  reduction  in  u>  leads  to  instability  of  the  numerical  method.  We  conclude 
that  the  amount  of  damping  is  a  strong  function  of  implicitness  parameter.  The 
results  shown  in  Fig.  2  were  for  a  sound  wave  with  a  cell  size  of  2.5  cm  using  a 
timestep  of  0.1  msec. 

The  dispersion  relation,  obtained  directly  from  the  calculations,  is  shown 
in  Fig.  3a  for  CFL  —  0.5.  The  CFL  number  is  defined  as 

sound  speed  X  time  step 

Ct  L  —  - - - . 

ceil  size 

These  calculations  were  made  by  varying  the  timestep  as  well  as  the  number  of 
cells  in  the  1  m  pipe.  The  product  of  the  wave  number,  k,  and  the  cell  size  Ax 
is  inversely  related  to  the  accuracy  of  representation  of  the  wave.  The  number 
oi  cells  per  wavelength  is  given  by 


N  = 


_2tt_ 
k  Ax 


On  the  vertical  axis,  we  show  and  ,  the  observed  and  theoretical  frequencies 
of  the  wave.  A  totally  dispersion  free  algorithm,  in  which  ouj  =  ut,  would 
yield  the  45  degree  line  shown.  Curves  for  different  values  of  the  implicitness 
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parameter  are  presented.  For  comparison  purposes,  the  results  for  the  explicit, 
predictor-corrector  JPBFCT  method  [4]  are  included.  JPBFCT  requires  two 
applications  of  the  FCT  algorithm  at  each  timestep.  This  two-step  process 
makes  JPBFCT  second  order  in  time. 

Figure  3b  indicates  the  change  in  amplitude  of  the  wave  in  one  period.  The 
amplification  is  always  less  than  unity,  indicating  that  the  wave  is  damped.  If 
the  amplification  were  greater  than  unity,  the  calculation  would  be  unstable.  As 
expected,  damping  increases  when  the  method  is  more  implicit.  Poorly  resolved 
wavelengths  are  damped,  even  by  the  fully  explicit  method.  It  should  be  noted 
that  BIC-FCT  with  w  —  0.55  performs  nearly  as  well  as  the  explicit  JBPFCT. 
Values  of  o>  nearer  0.5  brings  results  of  the  two  methods  even  closer. 

Figures  4a,  4b  and  5a,  5b  give  the  dispersion  relation  and  damping  for 
CFL  =  2  and  CFL  =  10  respectively.  Representation  of  the  sound  wave 
deteriorates  more  rapidly  as  resolution  is  lost  for  these  CFLs.  Curves  for  CFL  = 
10  are  shorter  than  those  for  lower  CFL  because  the  timestep  becomes  too  large 
to  resolve  the  oscillatory  nature  of  the  wave. 

This  example  points  out  the  need  for  caution  when  attempting  to  resolve 
sound  waves  at  high  CFLs.  Poorly  resolved  wavelengths  are  strongly  damped 
and  cannot  be  adequately  represented.  However,  if  only  long  wavelengths  ire 
of  interest,  BIC-FCT  can  provide  a  substantial  gain  over  an  explicit  method. 

A  Two-Dimensional  Problem 

When  BIC-FCT  is  applied  in  two  dimensions,  the  same  basic  three-step 
procedure  is  used.  In  addition,  we  use  time  splitting  in  the  two  spatial  di¬ 
mensions  to  implement  the  explicit  FCT  predictor  and  energy  corrector  step. 
However,  for  the  method  to  work,  the  elliptic  pressure  change  equation  must  be 
solved  in  two  dimensions. 

The  solution  of  the  elliptic  pressure  change  equation  is  a  substantial  part  of 
the  computional  effort  at  each  timestep.  In  one  dimension,  the  finite-difference 
form  of  the  pressure  difference  equation  can  be  solved  efficiently  in  0(N)  opera¬ 
tions,  where  N  is  the  number  of  grid  points,  using  standard  tridiagonal  methods 
(for  example,  see  Roache  [19]).  In  two  dimensions,  it  is  important  to  have  an 
efficient  elliptic  solver,  and  preferably  one  that  is  not  limited  to  specific  types 
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of  problems  with  specific  boundary  conditions.  In  the  calculation  presented  be¬ 
low  we  use  a  multigrid  method,  MGRID  [18],  which  is  very  fast  and  requires 
O(NlogN)  operations.  This  method  is  suitable  for  the  parallel  processing  in 
pipelined,  parallel,  and  vector  computers.  It  is  straightforward  to  use  any  other 
suitable  elliptic  solver. 

The  two-dimensional  Cartesian  test  problem  was  selected  to  demonstrate 
the  ability  of  BIC-FCT  to  treat  nearly  incompressible  swirling  flows.  A  potential 
vortex  with  a  central  core  was  used  as  an  initial  condition  in  a  square  10  m  x 
10  m  region.  The  initial  conditions  correspond  to  the  analytic  solution  of  a  line 
vortex  with  diffusion  which  is  of  the  form  [20] 


c  r 

^ tangential  =  ~ 


where  c  and  u  are  constants.  The  flow  very  rapidly  adjusts  to  the  presence 
of  the  walls,  but  this  does  not  affect  the  flow  close  to  the  vortex  center.  In 
this  test,  a  stretched  40  X  40  grid  was  used  with  the  smallest  cells  10  cm  in 
size  placed  at  the  center  of  the  vortex.  The  maximum  velocity,  at  the  start  of 
the  calculation,  was  30  m/s.  A  conservative  timestep  of  1  ms  was  used.  This 
should  be  contrasted  to  the  60  to  120  (is  timesteps  required  for  stability  in  a 
fully  explicit  method.  In  this  nearly  steady-state  problem,  the  effects  of  pressure 
fluctuations  are  expected  to  be  negligible.  Therfore,  we  could  use  u  =  1,  the 
fully  implicit  method. 

For  flow  visualization  purposes,  the  lower  half  of  the  fluid  has  been  marked 
and  appears  as  the  dark  area  in  Fig.  6.  In  the  absence  of  diffusion  processes, 
either  physical  or  numerical,  this  interface  remains  sharp  as  the  fluid  rotates  at  a 
constant  velocity.  Figures  7  and  8  show  the  position  of  the  interface  after  50  and 
200  timesteps  respectively.  The  interface  between  the  marked  and  unmarked 
fluid  is  no  longer  sharp,  due  to  numerical  diffusion.  The  interface  remains  fairly 
sharp  outside  the  core  region. 

The  velocity  decay  is  given  in  a  more  quantitative  manner  by  the  scatter 
plots  shown  in  the  next  set  of  figures.  Tangential  and  radial  components  of 
velocity  are  plotted  as  a  function  of  distance  from  the  vortex  center.  Crosses 
denote  the  velocity  actually  obtained  from  the  program  and  the  solid  line  pro¬ 
vides  a  least  squares  fit  of  the  data  to  the  analytic  solution  of  the  vortex  with 
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diffusion  [20].  The  initial  condition  is  shown  in  Figs.  9a  and  b.  Figures  10a, 
10b,  and  11a,  lib  show  the  velocity  after  50  and  200  timesteps,  respectively. 
The  peak  tangential  velocity  decreases  due  to  numerical  diffusion.  However, 
the  effective  diffusion  coefficient  is  not  a  constant  either  in  space  or  time  which 
leads  to  an  imperfect  fit  of  the  data  to  the  analytic  solution.  Scatter  in  the 
tangential  velocity  at  the  same  location  is  due  to  the  nonuniform  retardation 
caused  by  varying  amounts  of  numerical  diffusion.  Since  the  flow  is  essentially 
incompressible,  nonzero  radial  velocities  are  generated. 

We  now  examine  the  time  it  takes  to  do  one  computational  timestep.  Ta¬ 
ble  1  shows  a  timing  comparison  between  BIC-FCT  and  the  standard  module, 
JPBFCT,  very  similar  to  that  described  by  Boris  [4].  In  fact,  the  explicit  FCT 
predictor  in  BIC-FCT  is  similar  to  the  corrector  step  of  JPBFCT.  The  table 
shows  that  the  computational  time  required  per  timestep  compares  extremely 
favorably  to  that  for  the  explicit  method,  especially  at  the  larger  grid  sizes. 

IV.  SUMMARY  AND  DISCUSSION 

In  this  paper  we  have  described  the  Barely  Implicit  Correction  method, 
BIC,  for  calculating  subsonic  flows.  As  pointed  out  by  Casulli  and  Greenspan, 
only  the  pressure  and  velocity  terms  in  the  momentum  and  energy  equations 
respectively  have  to  be  treated  implicitly.  This  is  sufficient  to  remove  the  sound- 
speed  limit  on  the  timestep.  We  then  manipulated  the  equations  to  yield  a  single 
implicit  equation,  which  is  solved  for  a  correction  to  an  explicit  predictor  step. 
BIC  can  be  used  with  any  spatial  differencing  scheme.  BIC-FCT  provides  the 
accuracy  of  the  high-order  monotone  flux-corrected  transport  method  but  allows 
the  large  timesteps  possible  with  an  implicit  method. 

A  number  of  test  problems  showed  that  BIC-FCT  maintains  the  desirable 
high-order  monotone  characteristics  of  the  explicit  FCT  algorithm.  First,  we 
showed  that  it  could  propagate  a  contact  discontinuity  as  well  as  the  two-step 
JPBFCT.  We  also  presented  a  two-dimensional  example  of  a  swirling  flow. 

The  implicitness  parameter,  w,  plays  an  important  role  in  BIC-FCT  when¬ 
ever  sound  waves  and  pressure  oscillations  are  important  in  the  solution.  Damp- 
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ing  is  negligble  for  long  wavelenths  and  timesteps  when  w  =  0.5.  When  sound 
waves  are  not  important,  <jj  can  be  set  to  unity. 

The  major  gain  is  that  the  timestep  is  no  longer  restricted  by  the  sound 
speed.  This  improvement  is  achieved  at  little  or  no  additional  cost  per  timestep. 
The  cost  of  solving  the  elliptic  equation  is  recovered  by  the  elimination  of  the 
half-step  calculations  in  explicit  FCT. 

In  two-dimensional  problems,  an  efficient  method  of  solution  of  the  elliptic 
pressure  equation  is  essential.  The  multigrid  technique  MGRID  used  here  is 
among  the  fastest.  However,  the  application  of  this  technique  to  even  modestly 
complicated  geometries  is  not  straightforward.  Unstructured  multigrid  methods 
[21}  should  provide  the  necessary  flexibility. 

Equations  (9),  (10),  and  (16)  are  in  conservative  form.  Since  the  pres¬ 
sure  correction  only  appears  as  a  gradient  in  the  velocity  correction,  vorticity 
generation  and  transport  are  unchanged  by  BIC.  Thus  vorticity  stays  a  local, 
convected  quantity. 

In  summary,  BIC-FCT  has  opened  up  the  possibility  of  doing  very  accurate, 
very  slow  flow  calculations  in  which  compression  is  important.  Future  compu¬ 
tational  directions  include  extensions  to  finite  elements  [22]  and  to  addition  of 
other  physical  processes  such  as  gravity,  viscosity,  and  chemical  reactions  to 
simulate  premixed  flames,  diffusion  flames,  and  turbulent  jets. 
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Table  1. 

Timings*  per  Step  of  BIC-FCT  and  JPBFCT 


20  x  20 

40  X  40 

80  x  80 

BIC-FCT 

explicit 

6.8  ms 

17.0  ms 

54.1  ms 

elliptic 

3.8 

8.4 

22.5 

other 

2.7 

5.9 

17.1 

total 

13.3 

31.3 

93.7 

per  point 

33.3  fis 

19.6  /xs 

14.6  ns 

JPBFCT 

13.6  ms 

33.9  ms 

108.1  ms 

per  point 

34.0  fis 

21.2  fis 

16.9  /xs 
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Figure  1  —  Density  profiles  of  a  propagating  contact  discontinuity  at  50 

step  intervals. 
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(a)  dispersion  relation 

Figure  4  —  Dispersion  and  damping  of  sound  waves  by  BIC,  CFL  =  2.0. 
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(a)  dispersion  relation 


Figure  5  —  Dispersion  and  damping  of  sound  waves 
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Figure  10  (Continued)  —  Scatter  plot  of  velocity  in  vortex  flow,  50th  step. 
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(a)  radial  component  of  velocity 
Figure  1 1  —  Scatter  plot  of  velocity  in  vortex  flow,  200th  step. 
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